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It has been demonstrated that in the presence of weak collisions, described by the Lenard- 
Bemstein collision operator, the Landau-damped solutions become true eigenmodes of the 
system and constitute a complete set [C.-S. Ng, A. Bhattacharjee, and F. N. Skiff, Phys. 
Rev. Lett. 83, 1974 (1999), C.S. Ng, A. Bhattacharjee, and F. Skiff, Phys. Rev. Lett. 96, 
065002 (2004). We present numerical results from an Eulerian Vlasov code that incorpo- 
rates the Lenard-Bemstein collision operator [A. Lenard and I. B. Bernstein, Phys. Rev. 
112, 1456 (1958)]. The effect of collisions on the numerical recursion phenomenon seen 
in Vlasov codes is discussed. The code is benchmarked against exact linear eigenmode 
solutions in the presence of weak collisions, and a spectrum of Landau-damped solutions 
is determined within the limits of numerical resolution. Tests of the orthogonality and the 
completeness relation are presented. 

PACS numbers: 52.20Fs, 52.25.Dg, 52.27.Aj, 52.35.Fp, 52.65.Ff, 52.65.Vv 
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I. INTRODUCTION 


The problem of kinetic damping of plasma waves in a Maxwellian plasma is one of the most 
fundamental concepts in plasma physics. It is well known that Landau-damped solutions 1 are not 
true eigenmodes. The true eigenmodes for a collisionless plasma were obtained by Van Kampen 2 
and Case 3 . The Case- Van Kampen modes constitute a continuous spectrum, which correspond 
to the infinite number of degrees of freedom of the system. At this level of description, each 
particle is considered to move independent of any effects from surrounding particles, and there is 
no coupling. In most situations of physical interest where the initial conditions are smooth, a broad 
and continuous spectrum of Case-Van Kampen modes is excited. The Landau-damped solutions 
emerge, in the long-time limit, as remnants due to the interference of the continuous spectrum of 
singular Case- Van Kampen eigenmodes. 

Lenard and Bernstein (LB ) 4 reconsidered the kinetic problem using an operator of the Fokker- 
Planck type 5 . They obtained an exact analytic solution with a dispersion relation that formally 
yields the Landau root in the limit of zero collisions. However, they did not discuss the nature of 
the spectrum or address the issue of completeness of the eigenmodes in the presence of collisions. 
In retrospect, this appears a bit surprising because one of the important features of the Lenard- 
Bemstein collision operator, unlike the Bhatnagar-Gross-Krook (BGK) collision operator 6 , is that, 
in the limit of zero collision, it leads to a problem in singular perturbation theory in velocity space. 
While both BGK and LB operators produce the Landau solution in the limit of zero collision, the 
impact on the spectrum is profoundly different in the case of the LB operator, which arguably is 
more physical than the BGK operator. 

The most recent and direct impetus for theoretical studies on the nature of the kinetic spectrum 
in the presence of weak collisions have come from the remarkable experiments and analyses of 
Skiff and co-workers 7,8 , who used laser-induced fluorescence techniques to measure perturbed ion 
distribution functions in a stable plasma at unprecedented levels of accuracy (by two orders of 
magnitude compared with previous measurements). They offered the important insight that the 
linear eigenmode spectrum in the presence of weak collisions is intrinsically discrete. While the 
Vlasov description smoothens particle discreteness to form a continuous kinetic fluid in phase 
space and allows particles to move essentially independent of each other except through their role 
in collectively supporting a self-consistent electric field, binary collisions, embodied in the Fokker- 
Planck or the LB operator, introduce particle discreteness in the problem, leading to the emergence 
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of a discrete spectrum from the continuum modes. While the transformation of continuous spectra 
to discrete spectra behavior has been known to occur in strongly collisional systems that obey fluid 
equations (such as hydrodynamics or resistive MHD equations), the weakly collisional kinetic 
problem has not received the attention it deserves. Following the work of Skiff and co-workers, 
Ng, Bhattacharjee, and Skiff (NBS) 9 ’ 10 , demonstrated that in the presence of weak collisions, 
the singular Case- Van Kampen continuous spectrum is completely eliminated, and replaced by 
a discrete and smooth spectrum of eigenmodes which, furthermore, constitute a complete set. 
The Landau-damped solutions emerge as true eigenmodes of the weakly collisional theory in the 
limit of zero collision 9 . The demonstration of completeness of the discrete eigenmodes by NBS 10 
appears to resolve a controversy as to whether the total spectrum consists of a discrete as well as a 
continuous part 11 or only a discrete spectrum. It should be borne in mind, however, that the actual 
differentiation of the two types of spectra in real or numerical experiments is a subtle issue. 

The main goal of this paper is to test numerically the theoretical results obtained by NBS. For 
the first time, we report results from a kinetic Eulerian code in ld-lv space (that is, phase space 
consisting of one spatial coordinate x and one velocity space coordinate v) that includes the LB 
collision operator. We describe the code, hereafter referred to as the Kinetic Code, in Section II. 
In Section II E, we report on some standard tests of the code and compare the predictions of the 
code in the presence of collisions with the linear eigenmode analysis of NBS. In Section III, we 
discuss our effort in testing the completeness relation of NBS eigenmodes, which is a significant 
numerical challenge. We conclude in Section IV with a summary. 

II. KINETIC CODE 

A. Description of the Numerical Method 

The ld-lv Kinetic Code integrates the Vlasov-Poisson system, and has been extended to im- 
plement the Lenard-Bemstein collision operator in this work: 


Here / 


df_ 

dt 


d _l 

V dx 


qEdf d [ 2 df 

— JT = U JE v f + v thjr 
m ov ov ov 


V • E = — [ [f(x, v, t ) - fi(v, 0)] dv. 

^0 J —oo 


(1) 

( 2 ) 


f(x,v,t) is the electron distribution function, 0) is the spatially uniform ion 
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background, m is the electron mass, v t h is the electron thermal velocity, v is assumed to be a 
constant collision frequency, E is the self-consistent electric field and e 0 is the permittivity of free 
space. 

The system is typically initialized with a Maxwellian distribution, Eq.(3). At t — 0, this 
equilibrium is perturbed by a sinusoidal perturbation, Eq.(4) 




n 0 - 

Vzirvth 



f^\x,v) = esin (kx)f ( '°\x,v), 


(3) 

(4) 


which corresponds to a perturbed electric field, 

E^\x) = e^^cos (kx). (5) 

e 0 k 

Here e is a positive constant, measuring the size of the initial perturbation. 

The Kinetic Code employs a numerical method similar to the one proposed by Schumer and 
Holloway 12 , though it uses a finite difference scheme in space rather than a Fourier spectral 
method. The spatial boundaries are periodic. The code uses the split-step method of Schumer and 
Holloway with a 4th order Runge-Kutta method as the time integrator. For the velocity-dependent 
part of the distribution function we choose the symmetrically weighted Hermite functions as a 
spectral basis, defined as 

^n(v) = 0 n {v ) = C n e~ v 2 / 2 H n (v), ( 6 ) 

where H n (v) is the n th Hermite polynomial and C n = 7r -1 / 4 (2 n n!) -1 ' /2 is a normalization constant. 
These functions are orthonormal and satisfy the following recursion relations: 

vr (v) = 0^V +1 + (7) 

£rM = + 4 r-'M- oo 

Following Schumer and Holloway 12 , these recursion relations are used to numerically imple- 
ment the l.h.s. of Eq. (1). The spatial derivative is generally implemented using 2nd order central 
difference and the Poisson equation is solved by inverting a tri-diagonal matrix. However, for the 
simulations in this work, we used the linear, single-mode option of the code, which implements 
spatial derivatives exactly by multiplication with ik. 
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B. Implementation of the Lenard-Bernstein Collision Operator 


The Lenard-Bernstein collision operator 4 operator comes from the full Fokker-Planck opera- 
tor by linearizing and then assuming constant collision frequency, and single spatial and velocity 
dimensions. Its effect is to relax a non-thermal distribution function to a Maxwellian with a pre- 
scribed thermal velocity. It is known to conserve particle number and, for non-drifting distribution 
functions, momentum. These collisions have been implemented in the Kinetic Code by exploiting 
the recursion relations, Eqs. (7) and(8). 

Given a representation of the distribution function expanded in symmetrically weighted Her- 
mite functions, and leaving out the spatial and temporal dependence for brevity, 


/ 0 ) = ^2 fm'lpniv ), 


n = 0 


the right hand side of Eq. (1) is expanded and Eqs. (7), (8) are applied. 


(9) 
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( 10 ) 


Expanding the l.h.s. as well, 


LB = y^LBn^v), 

n = 0 


(ID 


and matching coefficients on both sides, we find the algorithm implemented in the code: 


LB n = v 


yj{n - 1 )n 


(1 - 4 ) 



2n + 1 
2 



+ !)(« + 2) (1 + <&) 


( 12 ) 
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FIG. 1 . Evolution of non-Maxwellian distribution function in the collisional Vlasov-Poisson system. Shown 
is the distribution function f(v) as a function of velocity at the initial time t = 0 (dashed), at an intermediate 
time t = 2 (dotted) where the oscillations have already strongly damped, and at t = 20 (solid), where the 
system has fully relaxed to a Maxwellian distribution. 

C. LB Operator Test: Thermalization of a Non-thermal Distribution 

The LB collision operator thermalizes non-Maxwellian distribution functions. Here we demon- 
strate how our implementation in the Kinetic Code exhibits this behavior. 

We start with a filamented perturbed Maxwellian initial condition for the distribution function 

f = f(v,t ): 

1 1,2 

f (u, 0) = (1 T 6 cos Kv) — -j = — e 2v th. (13) 

a/27 TVth 

Figure 1 shows the temporal evolution of the distribution function computed by the Kinetic 
Code, where we chose e = 0.3, K = 3, v t h = 3 and v = 0.01. It can be clearly seen how the 
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oscillatory perturbation gets rapidly damped and the system evolves to a Maxwellian distribution 
with width v t h- 

D. Effects of collisions on numerical recursion 

Eulerian Vlasov solvers are limited numerically by velocity-space filamentation of the distribu- 
tion function. In the absence of collisions, velocity shearing in phase space occurs inevitably due 
to the free-streaming or ballistic motion of the particles, represented by the first two terms in the 
left-hand- side ofEq. (I) 13-16 . 

This standard phenomenon is illustrated in Fig. 2. We begin at (a) t = 0 showing the initial 
perturbation in the distribution function j\ (see Eq. (4)). A little while later, in (b), particles with 
large velocities will have moved further in space than those with smaller velocities. At a later time 
in (c), the distribution function consists of fine filaments in velocity space. In the collisionless 
Vlasov-Poisson system, there is no mechanism to suppress the increasingly fine filamentation, so 
eventually the velocity scales fall below the grid size, leading to the well-known phenomenon of 
numerical recursion. In the context of the present method, one can postpone this time by increasing 
resolution or, to some extent, by rescaling the velocity space basis function, as discussed at length 
by Schumer and Holloway 12 . 

Figure 3 demonstrates the consequence of numerical recursion: The electric field in the col- 
lisionless run (thin line) follows the exponentially damped standing wave correctly initially, but 
recursion occurs at t ~ 150 where the simulated electric field shows erroneous behavior. The 
same plot demonstrates that a small collisionality of u = HP 4 is sufficient to avoid numerical 
recursion, so that the exponential decay continues. It should be noted that while the initial evo- 
lution is the same in both cases, in the collisional case the least damped collisional eigenmode 
eventually dominates as other modes in the initial condition decay away, which, by virtue of being 
an eigenmode, does not exhibit any further steepening of the velocity space gradients. Figure 4 
shows the structure of this eigenmode. 

To more closely investigate the occurence of recursion in Vlasov solvers and its suppression by 
the Lenard-Bernstein collision operator, it is instructive to consider the free- streaming approxima- 
tion to the Vlasov equation: 


df df d 
— — b r— = v— 


dt dx dv 



( 14 ) 
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FIG. 2. Evolution of the perturbed distribution function f\. As indicated in (a) by the arrows, different 
parts of the distribution function will move at their own velocity, leading to the distribution function being 
sheared and formation of filaments as shown at later times in (b) and (c). 


Setting the collision frequency u to zero for the moment, the equation is easily integrated along 
its characteristic, and in particular for a typical perturbation with wave number k in position space 
and Maxwellian distribution in velocity space f(x, v , 0) = a k cos (kx) exp(— v 2 / (2v^ h )), we find 
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FIG. 3. Evolution of electric field in the collisionless case (v = 0, light line), showing numerical recursion 
at t ~ 180 Wpg 1 and suppression of the recursion by collisionality u = 10~ 4 (thick line). 


the solution: 

f(x, v, t ) = f(x — vt, v, 0) = ak cos k(x — vt ) e 2v th . 


( 15 ) 


Taking a cut in velocity space at an arbitrary location xq, we see that the distribution function 
now has an oscillatory dependence oc cos (ktv) on v in addition the Maxwellian envelope, similar 
to our test initial condition in Fig. 1. The period in t'-spacc, 27 t/ (kt) decreases in time, and it is 
clear that the oscillation will become unresolved as we resolve a single cosine period with less 
than 4 points in v space: 


2vr 7T 

~ —— A:l\V T reC ur — " T . 

kT recur 2k Av 


( 16 ) 


As the width of Hermite functions is approximately the square root of their order, the resolution of 
our velocity-scaled basis functions is Av ~ , hence we obtain the recursion time as Schumer 
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FIG. 4. Structure of the late-time perturbed distribution function f\ in the collisional case, u = 10 -4 . The 
initial evolution is virtually identical the the collisionless evolution shown in Fig. 2, but at late time (shown 
is t = 253ce“g 1 ), the system has evolved into the least-damped collisional eigenmode. 


and Holloway 12 : 


TTy/N^ 

2kU ' 


(17) 


We performed a number of runs varying the velocity space resolution N u from 64 through 2048 
and determined the recursion times. Figure 5 confirms the expected Nu 2 scaling of the recursion 
time. 

While recursion at a finite time is unavoidable for numerical solutions of collisionless Vlasov- 
Poisson system, Ng et al. 9 showed that the character of the solutions changes fundamentally in 
the presence of collisions. In particular, the Landau damped solutions are transformed into actual 
eigenmodes, which as such remain unchanged during the temporal evolution, instead of becoming 
more and more oscillatory in velocity space. It is therefore expected that the numerical solver 
can follow the time evolution for long times, the phenomenon of recursion should disappear given 
sufficient resolution. 

Turning back to Eq.(14), the initial evolution of the system at finite but small collisionality fi 
will be virtually identical, since the initial Maxwellian profile does not have large gradients, so 
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FIG. 5. Numerical recursion time T recur as a function of velocity space resolution N u . The observed recur- 
sion times compare well to the model T reC ur ~ TTy/N u /(2kU). Parameters were k = 1 ,vth = 1/3, U = 
•5 v t h- 

the collisional terms on the right hand side are small. As time goes on, we again get the filamen- 
tation, leading to growing gradients in velocity space. Eventually, the collisional terms becomes 
important and serves to suppress growth of even higher modes, the dissipative term vv 2 h d 2 f/dv 2 
making the most important contribution. Substituting in our previously obtained free-streaming 
solution for / (Eq. (15)) and postulating that the collisional term quenches further filamentation 
once its contribution is comparable to the decay rate |t|, we find 

r “" “ t/?(£SF- <18) 

Recursion is prevented if further filamentation is suppressed before the recursion time, ie., 
Tcoii < T recur . To test this criterion, we solve for the minimum collision frequency needed to 
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FIG. 6. Collisionality v norecU r required to prevent recursion as a function of velocity space resolution N u . 
Parameters were k = 1 ,vth = 1/3, U = .5 vth- 


prevent recursion given a velocity space resolution N u : 

„ ( 2^) 2 7 


"norecur 


( 19 ) 


7T 2 (kv th ) 4 N u 

Figure 6 shows the collision frequency v needed to avoid the recursion previously observed in the 
collisionless system. It confirms the expected N~ l scaling. 


E. Kinetic Code Verification 

We have benchmarked the collisionless version of the Kinetic Code with roots of the plasma 
dispersion relation from a solver written by C.-S. Ng (referred to hereafter as the Plasma Disper- 
sion Function (PDF) Solver). The relevant dispersion equation is: 

l + a(l + fiZ(ft)) =0, (20) 
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FIG. 7. The Langmuir Dispersion Relation. The upper plots show the frequency of oscillation normalized 
by plasma frequency. The lower plot shows the damping rate normalized by plasma frequency. ’’PDF root” 
is the roots from the Plasma Dispersion Function Solver. ’’Code root” is the result from the Kinetic Code. 

where a = u 2 /{kv th ) 2 , C = u/ (\/2kv lh ) , and Z(Q) is the plasma dispersion function. As shown 
in Fig. 7, very good agreement is found between the Kinetic Code and the PDF Solver. 

We have also performed benchmarks of the damping of collisional Langmuir waves with the 
Kinetic Code. Table II, show that the damping rate matches the NBS eigensolver predictions. In 
the long-wavelength approximation, the eigenfrequencies, Q, take the form 

n = n 0 + n lf i, (2i) 


where, 


Ci = 


iVLq 


3 a V3 'at 


1 

20 ? 


-±+(l + ±) 

a \ a.) 




( 22 ) 
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TABLE I. NBS modesolver damping rates compared with the correction formulae, Eqs. (21) and (22). 


a 

fj, Dj NBS 

Uj Correction Formula 

Relative Error 

4.0 

0.01 -0.22254755 

-0.22253905 

0.96565 * 10 -5 

5.0 

0.01 -0.17044662 

-0.17044252 

0.97423 * 10" 5 

9.0 

0.1 -0.12701015 

-0.12863040 

0.96470 * 10 -3 

20.0 

0.1 -0.06674260 

-0.68091868 

0.55108 * 10" 3 

25.0 

0.1 -0.06149031 

-0.61999834 

0.35246 * 10" 3 


TABLE II. NBS Benchmark of Kinetic Code. Damping rates for various collision frequencies, // , recovered 
from fitting simulation results, compared to those predicted by the NBS eigenmode solver. 



n t NBS 

f Kinetic Code Relative Error 

0 

-0.0548864 

-0.0549166 

5.50 x 10” 4 

0.00001 

-0.0548937 

-0.0549240 

5.52 x 10 -4 

0.0001 

-0.0549601 

-0.0549904 

5.51 x 10 -4 

0.001 

-0.0556237 

-0.0556539 

5.43 x 10" 4 

0.01 

-0.0622458 

-0.0622758 

4.82 x 10" 4 

0.1 

-0.1270101 

-0.1269715 

3.04 x 10 -4 


Table I shows the comparison between the numerical results obtained by NBS, and the correction 
formulae (21) and (22). 


III. NBS EIGENMODE STUDY 

This section contains tests of the NBS eigenmodes. First, orthogonality of the modes is ex- 
amined analytically and numerically. Next, we show that the NBS eigenmodes behave in the 
collisional Vlasov-Poisson system as predicted. Finally, we present some evidence that are not 
consistent with the property of completeness. (For numerical reasons, discussed below, we are not 
able to test the completeness property definitively.) 
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A. Orthogonality 


NBS 9 renormalized and Fourier-transformed the distribution function / in space and time, and 
recast the collisional Vlasov-Poisson system, Eqs. (1), in terms of a normalized distribution func- 
tion g[y)\ 

( “ “ a)a ~ * /_„ Au ' )M = (“ s + It) ■ (23) 

where u = v/(\/2v 0 ), S'(w) = V2v 0 f/n 0 , g 0 = exp (—u 2 )/y/n, n 0 is the equilibrium electron 
density, g{u) = a/2dg 0 /du, and g = v/{g/2kv 0 ). This is an eigenmode equation with solutions 
g m belonging to complex eigenfrequencies Q m . NBS found the corresponding adjoint equation 
and showed that its adjoint solutions G m are related to the eigenfunctions g m by 

/»oo 

o (~V I 

G n (u ) cx g n (u)e u +—/= g n (u )du . (24) 

J — oo 

By construction, G„ and g m belonging to different eigenfrequencies f2 n , Q m are bi-orthogonal, 
that is, 


G n (u)g m (u)du = 0. 


(25) 


This condition allows us to determine the expansion coefficients c n for expanding an arbitrary 
function g(u) = J2m=n c mgm(u) in terms of eigenfunctions g m , after the G rn are appropriately 
normalized: 


Cm, 


G m {u)g[u)du. 


(26) 


We have confirmed numerically bi-orthogonality for a number of eigenfunctions obtained by 
the method described in 9 . The eigenfunctions are provided as a set of Hermite expansion coeffi- 
cients {a mn } such that 


OO 

g m (y) ^ ^ ^ m.n G n H ri ( V ) C 
n = 0 


(27) 


where C n = 1 / (tt 1//4 \/ 2 ;, 7c! ) . We express the adjoint function G rn in terms of Hermite polynomials, 

2 

too, but without the weight function e~ v to enable us to use the orthogonality property of Hermite 
polynomials to evaluate integrals. We write: 

OO 

G m {v ) = y ^A mn C n H n (v). (28) 

n = 0 
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TABLE III. Numerical confirmation of the orthonormality of the NBS eigenmodes. Shown is f G n g m du 
for Least Damped Mode and Modes 1 to 4 for g = 0.1 


n\m 

LDM 

1 

2 

3 

4 


LDM 12 3 4 

1 6.6788e-14 2.4654e-10 4.3470e-07 3.3390e-05 

1.0726e-14 1 1.9583e-09 2.5373e-06 1.6947e-04 

2. 1735e- 1 1 1.0749e-09 1 6.8720e-06 3.4759e-04 

2.667 le-08 9.6933e-07 4.7828e-06 1 7.6067e-04 

1.5850e-06 5.0094e-05 1.8717e-04 5.8853e-04 1 


Substituting the expansions for G n , g n into Eq. (24), multiplying by C m H m (u) exp(— u 2 ) and 
integrating over u in the usual manner, we determine the coefficients A nm : 


A 


nm 


bn [^nra "b & O'nO^mo] • 


(29) 


since g n (u)du = 7r 1//4 a n o- The b n are yet undetermined normalization constants. For any 
function g{u) expanded in Hermite functions, g(u) = a nC n H n (u) exp(— u 2 ), we can per- 

form a change of basis to find its coefficients {c m } for an expansion in collisional eigenfunctions: 


c m = G n (u)g(u)du, 


J — OO 
oo 


^ A-nk&k br, 


k = 0 


(1 + O') &n0 a 0 + E Q'nkQ'k 


k = 1 


(30) 

(31) 


The normalization coefficients can now be determined by plugging in g n for g, for which we 
require 


1 = 


C'nQndu, b r 


(1 + °0 a nO + Y. * 


l nk 


k= 1 


-1 


bn 


(1 + a ) a n0 + Y, 


2 


k = 1 


(32) 

(33) 


Table III shows an example calculation where we have calculated and normalized the adjoint 
functions G n and checked that the resulting functions are in fact bi-orthogonal. It is evident that 
the numerical eigenfunctions exhibit large and rapid oscillations near the phase speed of the wave, 
which makes it numerically challenging to determine them to very high accuracy. 
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FIG. 8. Time evolution of the k = 1 component of the electric field for separate simulations initialized with 
LDM and Modes 1 to 4, p = 0.1. Solid lines show the best fit to the exponential decay phase. 

B. Structure of the Eigenmodes 

The NBS eigensolver identified the least damped mode and 4 further modes for // = 0.1. The 
complex eigenfrequencies are given in Table IV. 

Figure 8 shows the time evolution of Vlasov simulations initialized with eigenmodes found 
by the NBS eigensolver for fi = 0.1. Indicated are the electric field (symbols) and least square 
fits to the exponential decay part of the evolution. It can be clearly seen that all modes initially 
show exponential decay as expected for an eigenmode, and examination of the full distribution 
function output confirmed this, as well. For all modes other than the least damped mode (LDM), 
the exponential decay eventually breaks down and we see further slow decay at the growth rate of 
the LDM. 

We compare the growth rates found from the Vlasov simulations to the imaginary part of Q 
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TABLE IV. Growth rates of the collisional eigenmodes at // = 0.1, showing that the results of the Vlasov 
code agree very well with the eigenvalues found by the NBS solver. 


Mode flj NBS Q, Vlasov code Relative Error 


LDM -0.12701014 -0.12700989 

1.97 x 10" 6 

1 -1.63253445 -1.63253647 

1.24 x 10" 6 

2 -2.44408929 -2.44409574 

2.64 x 10" 6 

3 -3.06241689 -3.06242860 

3.83 x 10" 6 

4 -3.58772051 -3.58797131 

6.99 x 10" 5 


calculated by the eigensolver in Table. IV. The Vlasov code clearly captures the temporal evolution 
of the eigenmodes to high accuracy in the initial phase of the simulations. 


The reason why the higher modes seemingly cease to evolve according to their linear dynamics 
can be understood better by considering the shape of the eigenfunctions in velocity space, as 
plotted in Fig. 9. All the eigenfunctions g m are normalized with respect to their zeroth moment, 
ie., / x ‘ x g m {v)dv = 1, hence they all contribute equally to the charge density and also the electric 
field. In particular, while all the modes in Fig. 8 start at the same value of the electric field, the 
distribution function g A of Mode 4 is actually up to 10 6 times larger than the distribution function 
of the least damped mode giDM- Due to numerical discretization error, the initial Mode 4 is not 
exact, but contains small contributions from other modes, in particular the least damped mode. An 
LDM error component of 10 -6 of the magnitude of Mode 4 would create an error in the electric 
field of order unity already. While the initial LDM error component is actually smaller, about 
10 -10 , it decays away much slower than the main mode, so eventually the electric field contributed 
by this small component overwhelms the electric field of the main mode and shows up as the 
slowly decaying evolution after the initial exponential phase. It should be noted that the structure 
of the eigenmodes is numerically rather challenging - they contain large oscillatory components 
that almost cancel out as one integrates over v- space, which leads to a significant loss of precision 
when performed in floating point arithmetic. 
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FIG. 9. Real part of the collisional eigenfunctions g m (v) shown on a semi-logarithmic scale. Solid lines 
indicate positive values, while dashed lines indicate negative values. Collisionality is // = 0.1. Note that 
the maximum amplitude of the modes grows with mode number. In the plot here, the (jldm has the lowest 
maximum amplitude and g 3 has the largest. 

C. Decomposition into Eigenmodes 

Our goal in this work is to analyze simulation results by decomposing the plasma evolution 
into a linear combination of eigenmodes, which all evolve according to their eigenfrequencies, 
explaining the dynamics of the system as interference of those eigenmodes. In order to test our 
method, we have performed a simulation that starts with a superposition of 5 modes, all them 
equally weighted with c m = 0.2. Figure 10 should be compared to Fig. 8 , only that this time we 
perform a single simulation evolving all 5 modes at once, and instead of focusing on the electric 
field, we decompose the distribution function back into its eigenmode basis coefficients c n (t). We 
reproduce fundamentally the same behavior, as one would expect for a linear system, but with 


20 



FIG. 10. Evolution of a superposition of LDM and Modes 1 to 4, initially weighted equally with coefficient 
0 . 2 . 


some caveats. LDM and modes 1 and 2 show the same evolution as before. Mode 3 turns over into 
the LDM behavior at an early point, while Mode 4 shows bleeding from other modes right away 
and is therefore not plotted. The main reason for the additional numerical difficulty lies in the fact 
that the numerical modes are not exactly orthogonal, as previously shown in Table III, which makes 
it impossible to obtain an exact decomposition. In addition, the eigenmodes and orthogonality 
relations are written for modes expanded of asymmetrically weighted Hermite functions, while the 
Kinetic Code uses symmetrically weighted polynomials, incurring additional conversions when 
setting up the initial conditions and analyzing distribution functions. In fact, the eigenmodes are so 
sensitive to small numerical errors that 64-bit floating point precision turned out to be insufficient 
to maintain even the already limited orthogonality under conversion between the bases. We have 
hence implemented the conversion routines in high-precision arithmetic using the mpmath Python 
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FIG. 11. Simulation initialized with a standing Langmuir wave, shown is the electric field E and the 
magnitude of the decomposition coefficients \c n \ vs time. 

module 17 . 

D. Collisional Landau damping of a standing wave 

We now analyze the temporal evolution of a standing Langmuir wave in the presence of col- 
lisions. The electric field, as shown in Fig. 11 shows the expected standing wave oscillation, 
enveloped by exponential decay at the expected growth rate from the least damped mode. The 
composition of the initial Maxwellian into LDM and Modes 1 to 4 is also shown, and those com- 
ponents behave similarly to the test case where we used an equally weighted superposition of 
modes as initial condition. Modes 1 to 4, as expected, show rapid exponentially decay, leaving 
only the least damped mode to support the observed total electric field. The evolution, in fact, 
looks quite similar to our synthetic test initial condition, where we superposed 5 modes with equal 
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weights, which subsequently decayed at their respective growth rates. In particular, the turning 
over of the fast exponential decay for the higher modes into slower decay at LDM rate is just 
the previously described numerical artifact. Our simulations also explain why the first “bump” 
in the electric field is higher than the exponential fit would have us expect, a phenomenon also 
observed in collisionless simulations, e.g. 12 . Landau damping describes the late time evolution 
of the electric field, however, there is an initial transient phase. In the context of the collisional 
system, the transient phase can be explained in terms of the additional modes which contribute 
to the initial condition. These modes also carry electric field, which explains the initially higher 
values, however, as they rapidly decay away, their contribution becomes insignificant quickly and 
the evolution of the electric field is well described by just the least damped mode itself. 


E. Completeness 


The NBS modes were shown to be a complete set of eigenmodes by 10 . First, we will examine 
the completeness problem by using numerically constructed linear eigenmodes from the method 
described at the end of the last section. Using this method, an eigenmode g m (v) is expanded in the 
form 

OO 

^ ^ ■ ( 34 ) 

n=0 

based on normalized variables. Numerically of course the summation over n can only be done up 
to a finite number of modes number, say n max , since only a finite number of a mn can be calculated. 
In the same way, the adjoint function corresponding to g m (v) can be obtained, 


G m {y') — ^ ] A mn H n (v). 


71=0 


The orthonormal condition requires 


9m{y')G n {v')d,V &mni 


(35) 


(36) 


or 

OO 

' &ml-A n l Smn- (37) 

z=o 

In other words, the transpose matrix A 1 is the right inverse of the matrix a. Numerically this 
relationship can be confirmed to a high accuracy by including enough modes in the summations 
in Eqs. (34) and (35), i.e., with a large enough n max . 
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Using the same formulation, the completeness condition then requires 


OO 


y, g m {v)G m {v') = 8{v- v'), 


( 38 ) 


m = 0 


or 



(39) 


In other words, the transpose matrix A 1 is also the left inverse of the matrix a, or that a is invert- 
ible. Again, numerically we can only sum up a finite number of terms in Eq. (39), or that m can 
only take up to m max , the number of eigenmodes that can be evaluated accurately. 

Unlike the orthonormal condition, Eq.(37), the completeness condition, Eq. (39), is much more 
difficult to test numerically. This is because the eigenmode g m (v) gets more and more singular in 
the v space for larger m, and thus is more and more difficult to be expanded accurately in the form 
of Eq. (34) for a given finite number of terms, n max . In other words, the number of eigenmodes, 
m max , that can be calculated accurately is not large enough. As a matter of fact, this is much 
more so for the temporal problem than the spatial problem. When m max is not too large, there are 
enough terms in the sum in Eq. (39) to get to converged values. 

Facing this great difficulty, one way to proceed in order to at least testing the numerical frame- 
work, rather than studying physically interesting cases, is to consider the spatial problem when 
the collision is strong. This is because g m (v) is more singular, i.e., with very fine structures in the 
v space with very small widths in the small /i limit. In the case of large g, the eigenmodes are 
smooth enough in the v space that a large m max of eigenmodes can be calculated accurately. 

As an example, we show in Table V, some results of the calculation of T pq as defined in Eq. (39) 
for the first few values of p and q, for the case with g — 1, a = 0.1, m max = 5020, n max = 
21600.We see that generally these values of T pq differ from the expected value of 5 pq in the order 
of O(10” 2 ).Note that Y pq is identically zero by symmetry if p and q are not both even or both 
odd. The reason we cannot get even better agreement is again because of a finite m max that can be 
calculated. In fact, in the calculation of each T pq , the trend in convergence is not consistent with 
the theoretical expectation, although the rate of convergences appears to be quite slow. 

IV. CONCLUSION 

In this paper, we report results from a ki netic Eulerian code in ld-lv space that includes the 
Lenard-Bernstein collision operator. Some standard tests of the code are given. We have shown 
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P q Tpg 


0 

0 

0.99999925 + 0.011) 

1 

1 

1.0000007- 7.1 x 10 -7 ) 

0 

2 

0.016 - 0.008) 

2 

0 

1.1 x 10" 6 - 0.008) 

2 

2 

0.989 - 0.0056z 

1 

3 

-3.4 x 10 -6 - 9.2 x 10 -7 ) 

3 

1 

-1.5 x 10 -6 + 2.6 x 10 -7 ) 

0 

4 

-0.028 + 0.0069) 

4 

0 

-2.6 x 10" 6 + 0.0069) 

3 

3 

1.000003 + 3.9 x 10~ 6 ) 

0 

6 

0.038 - 0.0063) 

6 

0 

5.0 x 10 -6 - 0.0063) 

2 

4 

0.0195 - 0.0049) 

4 

2 

0.0098 - 0.0049) 

1 

5 

6.5 x 10 -6 + 5.8 x 10" 6 ) 

5 

1 

2.3 x 10" 6 + 3.9 x 10" 7 ) 


TABLE V. Numerical Test of Completeness 

Numerical values of F pq as defined in Eq. (39) for the first few values of p and q, for the case with /j = 1, 

a = 0.1, m max = 5020, n max = 21600. 


that a non-thermal distribution is thermalized by the LB operator. We have examined the effect 
that collisions have on numerical recursion. Significantly, we have benchmarked the code with the 
NBS modes. 


We have discussed in detail our efforts in testing the orthogonality and completeness relation 
of NBS eigenmodes. The numerical challenges in providing definitive tests are quite formidable, 
but we have presented some evidence in support of these properties. 
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